library(ncdf4)
library(dplyr)
library(raster)
library(data.table)

#####################################################################################################
############################ 0. Prepare input data ##################################################
#####################################################################################################

#change to full local path to the replication data
wd1 <- "C:\\Users\\YFan\\OneDrive - NORCE\\NFood_replication\\"


#===============calc regional impacts===============
y.85 <- nc_open(paste0(wd1, "/data/yield_data/RCP85_crop_data.2006-2100.nc"))
region = c("NAM", "SEA", "EUA", "SAM", "CAM", "SAF", "Global")
box_n = c(60,   55,  60,   0,  30,  15, 90) #larger box for SEA and SAF
box_s = c(30,  -10,  36, -40,   0, -30, -90)
box_w = c(-130, 70, -10, -80,-110,  -12, -180) #do not change lon to 0:360, because it is difficult to clip SAF and EUA which cross lon=0 
box_e = c(-70, 130,  70, -35, -50,  50, 180) #re-organize the conditional matrix and sample matrix instead 

#create a conditional matrix to select grid samples for each region
v3 <- y.85$var[[3]] #read all info of a variable
v3$name;v3$units;v3$ndims #time is always the last dimension
varsize <- v3$varsize #dim sequence reversed for real data: lon=288 lat=192 cft=16 time=95
rm(v3) 

#===Important: after change lon and lonxy extent to -180:180, need to reformat gridmx so that grid sample locations in the conditional matrix are consistent with lonxy
#lon 0:178.75 shifted from position 1:144 to 145:288 in lon and from rows [1:144,] to [145:288,] in lonmx
lon <- ncvar_get(y.85,"lon")-180
lat <- ncvar_get(y.85,"lat")
latxy <- matrix(lat, nrow=length(lon), ncol=length(lat), byrow=T) #for matrix col is y, y is lat, construct latxy by lon
lonxy <- matrix(lon, nrow=length(lon), ncol=length(lat), byrow=F) #for matrix row is x, x is lon, construct lonxy by lat

gridmx <- array(1:(288*192),varsize[1:2]) #array of grid index/sample ids
#rotate grid id matrix, original rows [1:144,] corresponding to lonE 0:178.75 and [145:288,] corresponding to lonW -180:-1.25 should be swichted so that lon 0 locates at the center of map
gridmx.r <- rbind(gridmx[145:288,], gridmx[1:144,]) #now gridmx rows [1:144,] contain grid ids for lonW -180:-1.25, consistent with lonxy[1:144,1]

#When creating a raster from a matrix, the sense of rows and columns is opposite: raster sees row as y, col as x, but the matrix and netcdf array store lon(x) in row, lat(y) in col
#so need to transpose the matrix before plotting it as raster
plot(raster(t(gridmx.r), xmn=-180, xmx=180, ymn=-90, ymx=90))
plot(raster(t(latxy), xmn=-180, xmx=180, ymn=-90, ymx=90)) #upper-left corner x=0,y=192 [row=1,col=1] is for lat=-90, lon=-180  
plot(raster(t(lonxy), xmn=-180, xmx=180, ymn=-90, ymx=90)) #but maps and coordinates start from lower-left corner x=0,y=0 for lat=-90, lon=0/-180

#===NOTE:lonxy, latxy, gridmx, and netcdf data matrix all start from upper-left corner, i.e. x=0,y=192 [row=1,col=1] stands for lat=-90, lon=-180
#but the lon/lat coordinates of regional boxes or global map all start from lower-left corner, i.e. x=0,y=0 for lat=-90, lon=0 or -180!
#so the actual climate and yield data (in matrix or array) needs to be rearranged to start from lower-left corner, so they can be clipped by regional boxes
#this can be done by flipping the latxy, gridmx (which stores the sample ids of climate and yield data)

#check 1: when directly read netcdf file into raster format, original data matrix [1:288, 1:192] from netcdf is 
#seen as x=0:360 (1:288), y=-90:90 (1:192) and it renames x to cols and y to rows [1:192, 1:288]; 
ssp2.r <- raster(paste0(wd1, "/data/yield_data/RCP45_crop_data.2006-2100.nc"), var="CROP_AREA", level=1)
NAvalue(ssp2.r) <- max(values(ssp2.r)); 
str(ssp2.r)
ssp2.rr <- raster::rotate(ssp2.r) #rotate extent to -180:180
plot(ssp2.rr) #lower-left corner x=0,y=0 for lat=-90, lon=-180

#now it is not possible to change the climate and yield data, but the grid sample ids can be rearranged so that regional clipping is possible
#so have to flip the latxy and gridmx (upside down) and transpose the rows/cols to be consistent with the lon/lat coordinate system of regional boxes
plot(raster::flip(raster(t(gridmx.r), xmn=-180, xmx=180, ymn=-90, ymx=90), "y"))  #flip over the y direction
plot(raster::flip(raster(t(latxy), xmn=-180, xmx=180, ymn=-90, ymx=90), "y")) #now becomes lower-left corner [1,1] for lat=-90, lon=-180

#===Solution: need to rev(latxy) and flip(gridmx) so that lat=-90 lon=-180 start from lower-left corner, not upper-left corner
#new latxy: 
latxy.r <- matrix(rev(lat), nrow=length(lon), ncol=length(lat), byrow=T) 
#new gridmx: there is no function to flip a matrix, an alternative is to
#convert a list of (x=lon, y=lat, z=gridmx) to raster, then convert back to matrix: 
grid.s <- list(x=lon, y=lat, z=gridmx.r)
grid.r <- raster(grid.s)
plot(grid.r)
#convert raster back to matrix and transpose x/y
grid.mx <- t(as.matrix(grid.r))
#double check:  
plot(raster(t(grid.mx), xmn=-180, xmx=180, ymn=-90, ymx=90)) #grid.mx is flipped over gridmx.r to start from lower-left corner
plot(raster(t(latxy.r), xmn=-180, xmx=180, ymn=-90, ymx=90)) #now lower-left corner is lat=-90 lon=-180
#===now the conditional matrix generated by lonxy, latxy.r, and grid.mx can be clipped by regional coordinate boxes

save(list=c("region", "grid.mx", "latxy.r", "lonxy", "box_w", "box_e", "box_s", "box_n"), file=paste0(wd1,"/data/RData/Fig.Ext1.regional-gridmx-forclipping.Rdata"))
         


#******************************************************************************
#=============================Regional statistics===========================
#******************************************************************************

#calc regional and global modeled yield difference
load(paste0(wd1,"/data/RData/MLR-pergridregression.Rdata")) #read MLR regression coefficients
load(paste0(wd1,"/data/RData/MLR-pergrid-log-wtm-effect-Bootstrap-1000.Rdata")) #boostrap of original yearly data (not moving average)
load(paste0(wd1,"/data/RData/Fig.Ext1.regional-gridmx-forclipping.Rdata"))
#use all grids including crop failures in yield,yield0,yield1,yield2 to calc modeled change in global yield
dif.8crop <- readRDS(file = paste0(wd1,"/data/RData/allclimate-yield-changedata-MLR.Rds")) #all grid samples including those yield=0 (set to 5th percentile of each crop)
dif.8crop[,length(unique(Sample)), by=.(Scenario)] #same sample size
dif.8crop[Year==2090,sum(area.wt), by=.(Scenario,Year)]

##read original global landuse data for SSP2 and SSP5
#load data prepared by Data1.GrowingSeasonClimate-Yield.R
load(paste0(wd1,"/data/RData/GridSample-SSP2-SSP5-Data-final.Rdata"))
df.ssp2.all[Year==2090,sum(area),by=.(Scenario,Year)] #this is real yield and area data for RCP45, based on its own area mask (area>1000ha)
df.ssp5.all[Year==2090,sum(area),by=.(Scenario,Year)]
df.ssp5.all[Year==2090,summary(yield),by=.(Scenario,Year)]
df.ssp2.all[Year==2090,summary(yield),by=.(Scenario,Year)] 
#=do not apply yield mask, use all samples including crop failure (yield=0) to calc luc.eff and tot effect on global crop production
df.lu.ssp2 <- df.ssp2.all[Year>=2020 & Year !=2100]
df.lu.ssp5 <- df.ssp5.all[Year>=2020 & Year !=2100] #ssp5 has same number of samples as dif.8crop
length(unique(df.lu.ssp2$Sample)); length(unique(df.lu.ssp5$Sample)) #SSP2 has much more samples, more dispersed crop areas than SSP5



library(forecast)

#===========regional clipping and statistics of climate effects============
sample.all <- NULL
eff.crop.all <- eff.tot.all <- data.table()
for (n in 1:7){ #6 regions + global 
  #clip data to different regions (all grid samples for global)
  #each box yields a [1:288, 1:192] TRUE/FALSE conditional matrix, which is used to select data
  print(paste("calculate for region:", region[n]))
  
  sample.reg <- grid.mx[latxy.r>=box_s[n] & latxy.r<=box_n[n] & 
                        lonxy>=box_w[n] & lonxy<=box_e[n]] 
  sample.all <- c(sample.all, list(sample.reg))
    
  #if (n < 7) {  
    
    #======option2: use grid-based regression coeffs to predict partial and total climate effects and then aggregate to each region
    #====estimate regional climate, co2, luc effects (log) 
    coeff.i.reg <- coeff.i.samples[Sample%in%sample.reg]
    coeff.r.reg <- coeff.r.samples[Sample%in%sample.reg] #grid sample specific coeffs

    effect.i.reg <- coeff.i.reg[dif.8crop.s[Irrig==TRUE & Sample%in%sample.reg],                                    
                                .(Intercept, T.log=T*i.tsa.dif, RD.log=RD*i.radd.dif, RI.log=RI*i.radi.dif, 
                                  R.log=RD*i.radd.dif+RI*i.radi.dif, P.log=0, RH.log=0, M.log=0, 
                                  co2.log=i.co2.log, luc.log=i.luc.log, 
                                  Y.mod.dif=i.yield.dif, Y.mod.dif2=i.yield.dif+i.co2.log, Y.mod.dif3=i.yield.dif+i.co2.log+i.luc.log, 
                                  yield1=i.yield1, yield0=i.yield0, yield=i.yield, 
                                  area=i.area.wt, Year=i.Year), by=.EACHI, on=.(Scenario,Crop,Sample)]
    effect.r.reg <- coeff.r.reg[dif.8crop.s[Irrig==FALSE & Sample%in%sample.reg], 
                                .(Intercept, T.log=T*i.tsa.dif, R.log=RD*i.radd.dif+RI*i.radi.dif, M.log=P*i.ppt.dif+RH*i.rh.dif,
                                  RD.log=RD*i.radd.dif, RI.log=RI*i.radi.dif, P.log=P*i.ppt.dif, RH.log=RH*i.rh.dif, 
                                  co2.log=i.co2.log, luc.log=i.luc.log,
                                  Y.mod.dif=i.yield.dif, Y.mod.dif2=i.yield.dif+i.co2.log, Y.mod.dif3=i.yield.dif+i.co2.log+i.luc.log, 
                                  yield1=i.yield1, yield0=i.yield0, yield=i.yield, 
                                  area=i.area.wt, Year=i.Year), by=.EACHI, on=.(Scenario,Crop,Sample)]
    effect.i.reg  <- effect.i.reg[,Tot.log:=T.log+RD.log+RI.log+Intercept] #if regression includes intercept, do not forget it here.
    effect.r.reg  <- effect.r.reg[,Tot.log:=T.log+RD.log+RI.log+P.log+RH.log+Intercept]
    effect.i.reg  <- effect.i.reg[,Tot.log2:=Tot.log+co2.log]  #add co2 effect
    effect.r.reg  <- effect.r.reg[,Tot.log2:=Tot.log+co2.log]
    effect.i.reg  <- effect.i.reg[,Tot.log3:=Tot.log2+luc.log] #all effects
    effect.r.reg  <- effect.r.reg[,Tot.log3:=Tot.log2+luc.log] #larger than adding global average luc.log to tot.log2 after aggregation
    effect.i.reg  <- effect.i.reg[,Res.log:=Y.mod.dif-Tot.log] #residuals
    effect.r.reg  <- effect.r.reg[,Res.log:=Y.mod.dif-Tot.log] 
    
    #merge rainfed and irrigated
    effect.i.reg$Irrig=TRUE
    effect.r.reg$Irrig=FALSE
    effect.reg <- rbind(effect.i.reg, effect.r.reg)

    #from grid sample to global average per crop
    #final option: weight average log effect by crop area (same as Bootstrap method)
    eff.reg.log <- effect.reg[, .(T.log=weighted.mean(T.log, area, na.rm = T), RD.log=weighted.mean(RD.log, area, na.rm = T),RI.log=weighted.mean(RI.log, area, na.rm = T),
                               P.log=weighted.mean(P.log, area, na.rm = T),RH.log=weighted.mean(RH.log, area, na.rm = T),
                               R.log=weighted.mean(R.log, area, na.rm = T), M.log=weighted.mean(M.log, area, na.rm = T),
                               Tot.log=weighted.mean(Tot.log, area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, area, na.rm = T), 
                               co2.log=weighted.mean(co2.log, area, na.rm = T),luc.log=weighted.mean(luc.log, area, na.rm = T), 
                               Tot.log2=weighted.mean(Tot.log2, area, na.rm = T), Y.mod.dif2=weighted.mean(Y.mod.dif2, area, na.rm = T),
                               Tot.log3=weighted.mean(Tot.log3, area, na.rm = T), Y.mod.dif3=weighted.mean(Y.mod.dif3, area, na.rm = T),
                               Res.log=weighted.mean(Res.log, area, na.rm = T),
                               yield1=weighted.mean(yield1, area, na.rm = T), yield0=weighted.mean(yield0, area, na.rm = T), yield=weighted.mean(yield, area, na.rm = T),
                               area=sum(area, na.rm=T)),
                           by =.(Scenario,Crop,Irrig,Year)]
    
    #25 year moving average to get mean effect of 2075-2099 (at year 2087)
    # eff.reg.log.ma <- eff.reg.log[, lapply(.SD, ma, order=25), by =.(Scenario,Crop,Irrig)] # year 2087 corresponds to mean of 2075-2099
    # eff.reg.log.ma <- eff.reg.log.ma[, lapply(.SD, as.numeric), by =.(Scenario,Crop,Irrig)]
    # eff.reg.log.ma[,Year:=as.integer(Year)]
    # eff.reg.log.ma <- eff.reg.log.ma[!is.na(Year)]
    eff.reg.log.25yr <- eff.reg.log[Year>=2075 & Year<=2099, lapply(.SD, mean, na.rm=T), by =.(Scenario,Crop,Irrig)]

    #===get reference data
    #regional area weighted mean of grid-level change (CLM5 modeled) from full dataset dif.8crop
    dif.8crop.reg <- rbind(dif.8crop[Sample%in%sample.reg & Scenario=="RCP45 - RCP85" ,
                                    .(Y.mod.dif=weighted.mean(yield.dif, area.wt, na.rm = T),
                                      Y.mod.dif2=weighted.mean(log(yield)-log(yield0), area.wt, na.rm = T), #==yield.dif+co2.log,
                                      Y.mod.dif3=weighted.mean(yield.dif+co2.log+luc.log, area.wt, na.rm = T),
                                      yield1=weighted.mean(yield1, area.wt, na.rm=T), yield0=weighted.mean(yield0, area.wt, na.rm=T), yield=weighted.mean(yield, area.wt, na.rm=T),
                                      co2.log=weighted.mean(co2.log, area.wt, na.rm=T), luc.log=weighted.mean(luc.log, area.wt, na.rm = T),
                                      area=sum(area.wt, na.rm=T)), 
                                    by=.(Scenario,Crop,Irrig,Year)],
                          dif.8crop[Scenario!="RCP45 - RCP85" & Sample%in%sample.reg,
                                    .(Y.mod.dif=weighted.mean(yield.dif, area.wt, na.rm = T),
                                      yield0=weighted.mean(yield0, area.wt, na.rm=T), yield=weighted.mean(yield, area.wt, na.rm=T), luc.log=0, co2.log=0,
                                      area=sum(area.wt, na.rm=T)),
                                    by=.(Scenario,Crop,Irrig,Year) ], fill=TRUE)
    
    #replace the crop area and base yield data of dif.8crop.s with global data for calculating Y.mod.dif3 (pct change of global average yield)
    lu.ssp2.reg <- df.lu.ssp2[Sample%in%sample.reg,  #the sample size here is different from dif.8crop.reg due to LUC
                                 .(yield2=weighted.mean(yield, area, na.rm=T), fertn2=weighted.mean(fertn, area,na.rm=T), 
                                   area2=sum(area, na.rm=T)), by=.(Scenario,Crop,Irrig,Year)]
    lu.ssp2.reg[Scenario=="SSP2",Scenario:="RCP45 - RCP85"]
    dif.8crop.reg <-dif.8crop.reg[lu.ssp2.reg, yield2:=i.yield2, on=.(Scenario,Crop,Irrig,Year)] #yield of RCP45(SSP2)
    dif.8crop.reg <-dif.8crop.reg[lu.ssp2.reg, area2:=i.area2, on=.(Scenario,Crop,Irrig,Year)] #area of RCP45(SSP2)
    
    #=not all regions have every crop type across years
    dif.8crop.reg.25yr <- dif.8crop.reg[Year>=2075 & Year<=2099, lapply(.SD, mean, na.rm=T), by =.(Scenario,Crop,Irrig)]
    
    #convert log to % effect (always after aggregation)
    #averaged trends
    # eff.reg.ma <- eff.reg.log.ma[dif.8crop.reg.ma,
    eff.reg.25yr <- eff.reg.log.25yr[dif.8crop.reg.25yr,
                                 .(T.eff=(exp(T.log) - 1)*100, R.eff=(exp(R.log) - 1)*100, M.eff=(exp(M.log) - 1)*100, #must merge RD+RI, P+RH in log term before aggregation
                                   RD.eff=(exp(RD.log)-1)*100, RI.eff=(exp(RI.log) - 1)*100, P.eff=(exp(P.log) - 1)*100, RH.eff=(exp(RH.log) - 1)*100,
                                   co2.eff=(exp(i.co2.log)-1)*100, luc.eff=(exp(i.luc.log)-1)*100, #luc.eff=(exp(luc.log)-1)*100, #LUC from dif.8crop (15.4%) is slightly higher than from dif.8crop.s (14.9%)
                                   Tot.eff=(exp(Tot.log) - 1)*100, Tot.eff2=(exp(Tot.log2) - 1)*100, #co2.eff=(exp(co2.log)-1)*100, #co2 effect from dif.8crop (-19.2%) is slightly higher than from dif.8crop.s (-18.3%)
                                   Tot.eff3=(exp(Tot.log3) - 1)*100, #climate+co2+LUC added at each grid sample in dif.8crop.s and then aggregated to crop (8.009%) is larger than adding global co2.log and luc.log after aggregation (7.89% or 6.96%)
                                   Y.mod.dif=(exp(i.Y.mod.dif) - 1)*100, Y.mod.dif2=(exp(i.Y.mod.dif2) - 1)*100,
                                   Y.mod.dif3=(exp(i.Y.mod.dif3)-1)*100, #wtm of per grid modeled change from dif.8crop (climate+co2+LUC added at each grid sample before aggregation) 7.13%
                                   #Y.mod.dif3=(exp(i.Y.mod.dif+i.co2.log+i.luc.log)-1)*100, #adding global co2.log and luc.log after aggregation (7.01%)
                                   #Tot.eff3=(exp(Tot.log+co2.log+luc.log)-1)*100, #adding co2.log and luc.log from dif.8crop.s gives slightly higher tot effect (7.89%) than from dif.8crop (6.96%)
                                   #Tot.eff3=(exp(Tot.log+i.co2.log+i.luc.log)-1)*100, #adding co2.log and luc.log from dif.8crop after aggregation to crop matches better with modeled global average change (6.96 vs. 7.01%)
                                   yield=i.yield,yield0=i.yield0,yield1=i.yield1,yield2=i.yield2, area=i.area, area2=i.area2
                                 ), by=.EACHI, on=.(Scenario,Crop,Irrig)]
    eff.reg.25yr[Scenario!="RCP45 - RCP85", co2.eff:=NA]
    eff.reg.25yr[Scenario!="RCP45 - RCP85", luc.eff:=NA]
    
   rm(effect.i.reg, effect.r.reg,effect.reg,dif.8crop.reg)
   

  #====merge effects of 8 crops to 6 crops and merge irrig/rainfed
  eff.pct.all<- copy(eff.reg.25yr)
  eff.pct.all[Crop=="Tropical Corn", Crop:="Corn"]
  eff.pct.all[Crop=="Tropical Soybean", Crop:="Soybean"]
  
  #per crop effect: merge irrig/rainfed, weight by area and base yield
  eff.crop <- rbind(eff.pct.all[Scenario=="RCP45 - RCP85", 
                               .(T.eff=weighted.mean(T.eff, yield1*area, na.rm = T), R.eff=weighted.mean(R.eff, yield1*area, na.rm = T),M.eff=weighted.mean(M.eff, yield1*area, na.rm = T),
                                 RD.eff=weighted.mean(RD.eff, yield1*area, na.rm = T),P.eff=weighted.mean(P.eff, yield1*area, na.rm = T),
                                 RI.eff=weighted.mean(RI.eff, yield1*area, na.rm = T),RH.eff=weighted.mean(RH.eff, yield1*area, na.rm = T),
                                 Tot.eff=weighted.mean(Tot.eff, yield1*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield1*area, na.rm = T),
                                 Tot.eff2=weighted.mean(Tot.eff2, yield0*area, na.rm = T), Y.mod.dif2=weighted.mean(Y.mod.dif2, yield0*area, na.rm = T), #climate+co2 effect
                                 Y.mod.dif3=weighted.mean(Y.mod.dif3, yield0*area, na.rm = T), #climate+co2+LUC
                                 Tot.eff3=weighted.mean(Tot.eff3, yield0*area, na.rm = T), #climate+co2+LUC added at each grid sample before aggregation
                                 co2.eff=weighted.mean(co2.eff, yield0*area, na.rm = T), luc.eff=weighted.mean(luc.eff, yield*area, na.rm = T), #luc.log = log(i.yield2)-log(i.yield)
                                 yield0=weighted.mean(yield0, area, na.rm=T),yield1=weighted.mean(yield1, area, na.rm=T),
                                 yield=weighted.mean(yield, area, na.rm=T), yield2=weighted.mean(yield2, area2, na.rm=T),
                                 area=sum(area, na.rm = T), area2=sum(area2, na.rm = T)), by =.(Scenario,Crop)],
                   eff.pct.all[Scenario!="RCP45 - RCP85", 
                               .(T.eff=weighted.mean(T.eff, yield0*area, na.rm = T), R.eff=weighted.mean(R.eff, yield0*area, na.rm = T),M.eff=weighted.mean(M.eff, yield0*area, na.rm = T),
                                 Tot.eff=weighted.mean(Tot.eff, yield0*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield0*area, na.rm = T), 
                                 RD.eff=weighted.mean(RD.eff, yield0*area, na.rm = T),P.eff=weighted.mean(P.eff, yield0*area, na.rm = T),
                                 RI.eff=weighted.mean(RI.eff, yield0*area, na.rm = T),RH.eff=weighted.mean(RH.eff, yield0*area, na.rm = T),
                                 yield=weighted.mean(yield, area, na.rm=T), yield0=weighted.mean(yield0, area, na.rm=T),
                                 area=sum(area, na.rm = T)), by =.(Scenario,Crop)], fill=TRUE)
  
  #compare with land use effect and total modeled change based on global average yield data
  # eff.crop <- eff.crop[Scenario=="RCP45 - RCP85", co2.eff:=(yield1-yield0)/yield0*100, by=.(Scenario,Crop,Year)] 
  # eff.crop <- eff.crop[Scenario=="RCP45 - RCP85", luc.eff:=(yield2-yield)/yield*100, by=.(Scenario,Crop,Year)] #aggregated differently from luc.dif2
  # eff.crop <- eff.crop[Scenario=="RCP45 - RCP85", Y.mod.dif:=(yield2-yield0)/yield0*100, by=.(Scenario,Crop,Year)] #aggregated differently from Y.mod.dif2
  # eff.crop <- eff.crop[Scenario!="RCP45 - RCP85", Y.mod.dif:=(yield-yield0)/yield0*100, by=.(Scenario,Crop,Year)] #same as weighting % effect by area*yield0 

  #get global total impact
  eff.tot <- rbind(eff.pct.all[Scenario=="RCP45 - RCP85", 
                               .(T.eff=weighted.mean(T.eff, yield1*area, na.rm = T), R.eff=weighted.mean(R.eff, yield1*area, na.rm = T),M.eff=weighted.mean(M.eff, yield1*area, na.rm = T),
                                 RD.eff=weighted.mean(RD.eff, yield1*area, na.rm = T),P.eff=weighted.mean(P.eff, yield1*area, na.rm = T),
                                 RI.eff=weighted.mean(RI.eff, yield1*area, na.rm = T),RH.eff=weighted.mean(RH.eff, yield1*area, na.rm = T),
                                 Tot.eff=weighted.mean(Tot.eff, yield1*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield1*area, na.rm = T),
                                 Tot.eff2=weighted.mean(Tot.eff2, yield0*area, na.rm = T), Y.mod.dif2=weighted.mean(Y.mod.dif2, yield0*area, na.rm = T), #climate+co2 effect
                                 Tot.eff3=weighted.mean(Tot.eff3, yield0*area, na.rm = T), Y.mod.dif3=weighted.mean(Y.mod.dif3, yield0*area, na.rm = T), #climate+co2+LUC
                                 co2.eff=weighted.mean(co2.eff, yield0*area, na.rm = T), luc.eff=weighted.mean(luc.eff, yield*area, na.rm = T), #luc.log = log(i.yield2)-log(i.yield)
                                 yield0=weighted.mean(yield0, area, na.rm=T),yield1=weighted.mean(yield1, area, na.rm=T),
                                 yield=weighted.mean(yield, area, na.rm=T), yield2=weighted.mean(yield2, area2, na.rm=T),
                                 area=sum(area, na.rm = T), area2=sum(area2, na.rm = T)), by =.(Scenario)],
                   eff.pct.all[Scenario!="RCP45 - RCP85", 
                               .(T.eff=weighted.mean(T.eff, yield0*area, na.rm = T), R.eff=weighted.mean(R.eff, yield0*area, na.rm = T),M.eff=weighted.mean(M.eff, yield0*area, na.rm = T),
                                 Tot.eff=weighted.mean(Tot.eff, yield0*area, na.rm = T), Y.mod.dif=weighted.mean(Y.mod.dif, yield0*area, na.rm = T), 
                                 RD.eff=weighted.mean(RD.eff, yield0*area, na.rm = T),P.eff=weighted.mean(P.eff, yield0*area, na.rm = T),
                                 RI.eff=weighted.mean(RI.eff, yield0*area, na.rm = T),RH.eff=weighted.mean(RH.eff, yield0*area, na.rm = T),
                                 yield=weighted.mean(yield, area, na.rm=T), yield0=weighted.mean(yield0, area, na.rm=T),
                                 area=sum(area, na.rm = T)), by =.(Scenario)], fill=TRUE)  
  # compare with land use effect and total modeled change based on global average yield data
  # eff.tot <- eff.tot[Scenario=="RCP45 - RCP85", co2.eff:=(yield1-yield0)/yield0*100, by=.(Scenario,Year)]
  # eff.tot <- eff.tot[Scenario=="RCP45 - RCP85", luc.eff:=(yield2-yield)/yield*100, by=.(Scenario,Year)]
  # eff.tot <- eff.tot[Scenario=="RCP45 - RCP85", Y.mod.dif:=(yield2-yield0)/yield0*100, by=.(Scenario,Year)]
  # eff.tot <- eff.tot[Scenario!="RCP45 - RCP85", Y.mod.dif:=(yield-yield0)/yield0*100, by=.(Scenario,Year)] #same as weighting % effect by area*yield0 
  
  eff.crop$Region=region[n]
  eff.tot$Region=region[n]
  eff.crop.all <- rbind(eff.crop.all, eff.crop)
  eff.tot.all <- rbind(eff.tot.all, eff.tot)
}

save(list = c("eff.crop.all","eff.tot.all"), file = paste0(wd1,"/data/RData/Fig.Ext1-2.Radar-regional-effect.RData")) #statistical prediction matches well with modeled difference


